home *** CD-ROM | disk | FTP | other *** search
/ STraTOS 1997 April & May / STraTOS 1 - 1997 April & May.iso / CD01 / INTERNET / SITES / LITTLE / P3SRC.ZIP / ATARI / QUADRICS.C < prev    next >
Encoding:
C/C++ Source or Header  |  1996-06-26  |  29.0 KB  |  1,543 lines
  1. /****************************************************************************
  2. *                quadrics.c
  3. *
  4. *  This module implements the code for the quadric shape primitive.
  5. *
  6. *  from Persistence of Vision(tm) Ray Tracer
  7. *  Copyright 1996 Persistence of Vision Team
  8. *---------------------------------------------------------------------------
  9. *  NOTICE: This source code file is provided so that users may experiment
  10. *  with enhancements to POV-Ray and to port the software to platforms other
  11. *  than those supported by the POV-Ray Team.  There are strict rules under
  12. *  which you are permitted to use this file.  The rules are in the file
  13. *  named POVLEGAL.DOC which should be distributed with this file. If
  14. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  15. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  16. *  Forum.  The latest version of POV-Ray may be found there as well.
  17. *
  18. * This program is based on the popular DKB raytracer version 2.12.
  19. * DKBTrace was originally written by David K. Buck.
  20. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  21. *
  22. *****************************************************************************/
  23.  
  24. #include "frame.h"
  25. #include "povray.h"
  26. #include "vector.h"
  27. #include "povproto.h"
  28. #include "bbox.h"
  29. #include "objects.h"
  30. #include "matrices.h"
  31. #include "planes.h"
  32. #include "quadrics.h"
  33.  
  34.  
  35.  
  36. /*****************************************************************************
  37. * Local preprocessor defines
  38. ******************************************************************************/
  39.  
  40. /* The following defines make typing much easier! [DB 7/94] */
  41.  
  42. #define Xd Ray->Direction[X]
  43. #define Yd Ray->Direction[Y]
  44. #define Zd Ray->Direction[Z]
  45.  
  46. #define Xo Ray->Initial[X]
  47. #define Yo Ray->Initial[Y]
  48. #define Zo Ray->Initial[Z]
  49.  
  50. #define QA Quadric->Square_Terms[X]
  51. #define QE Quadric->Square_Terms[Y]
  52. #define QH Quadric->Square_Terms[Z]
  53.  
  54. #define QB Quadric->Mixed_Terms[X]
  55. #define QC Quadric->Mixed_Terms[Y]
  56. #define QF Quadric->Mixed_Terms[Z]
  57.  
  58. #define QD Quadric->Terms[X]
  59. #define QG Quadric->Terms[Y]
  60. #define QI Quadric->Terms[Z]
  61.  
  62. #define QJ Quadric->Constant
  63.  
  64.  
  65.  
  66. /*****************************************************************************
  67. * Static functions
  68. ******************************************************************************/
  69.  
  70. static int Intersect_Quadric PARAMS((RAY *Ray, QUADRIC *Quadric, DBL *Depth1, DBL *Depth2));
  71. static void Quadric_To_Matrix PARAMS((QUADRIC *Quadric, MATRIX Matrix));
  72. static void Matrix_To_Quadric PARAMS((MATRIX Matrix, QUADRIC *Quadric));
  73. static int All_Quadric_Intersections PARAMS((OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack));
  74. static int Inside_Quadric PARAMS((VECTOR IPoint, OBJECT *Object));
  75. static void Quadric_Normal PARAMS((VECTOR Result, OBJECT *Object, INTERSECTION *Inter));
  76. static void *Copy_Quadric PARAMS((OBJECT *Object));
  77. static void Translate_Quadric PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  78. static void Rotate_Quadric PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  79. static void Scale_Quadric PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  80. static void Transform_Quadric PARAMS((OBJECT *Object, TRANSFORM *Trans));
  81. static void Invert_Quadric PARAMS((OBJECT *Object));
  82. static void Destroy_Quadric PARAMS((OBJECT *Object));
  83.  
  84. /*****************************************************************************
  85. * Local variables
  86. ******************************************************************************/
  87.  
  88. METHODS Quadric_Methods =
  89. {
  90.   All_Quadric_Intersections,
  91.   Inside_Quadric, Quadric_Normal,
  92.   Copy_Quadric,
  93.   Translate_Quadric, Rotate_Quadric,
  94.   Scale_Quadric, Transform_Quadric, Invert_Quadric,
  95.   Destroy_Quadric
  96. };
  97.  
  98.  
  99.  
  100.  
  101. /*****************************************************************************
  102. *
  103. * FUNCTION
  104. *
  105. *   All_Quadric_Intersections
  106. *
  107. * INPUT
  108. *   
  109. * OUTPUT
  110. *   
  111. * RETURNS
  112. *   
  113. * AUTHOR
  114. *
  115. *   POV-Ray Team
  116. *   
  117. * DESCRIPTION
  118. *
  119. *   -
  120. *
  121. * CHANGES
  122. *
  123. *   -
  124. *
  125. ******************************************************************************/
  126.  
  127. static int All_Quadric_Intersections(Object, Ray, Depth_Stack)
  128. OBJECT *Object;
  129. RAY *Ray;
  130. ISTACK *Depth_Stack;
  131. {
  132.   DBL Depth1, Depth2;
  133.   VECTOR IPoint;
  134.   register int Intersection_Found;
  135.  
  136.   Intersection_Found = FALSE;
  137.  
  138.   if (Intersect_Quadric(Ray, (QUADRIC *)Object, &Depth1, &Depth2))
  139.   {
  140.     if ((Depth1 > Small_Tolerance) && (Depth1 < Max_Distance))
  141.     {
  142.       VEvaluateRay(IPoint, Ray->Initial, Depth1, Ray->Direction);
  143.  
  144.       if (Point_In_Clip(IPoint, Object->Clip))
  145.       {
  146.         push_entry(Depth1, IPoint, Object, Depth_Stack);
  147.  
  148.         Intersection_Found = TRUE;
  149.       }
  150.     }
  151.  
  152.     if ((Depth2 > Small_Tolerance) && (Depth2 < Max_Distance))
  153.     {
  154.       VEvaluateRay(IPoint, Ray->Initial, Depth2, Ray->Direction);
  155.  
  156.       if (Point_In_Clip(IPoint, Object->Clip))
  157.       {
  158.         push_entry(Depth2, IPoint, Object, Depth_Stack);
  159.  
  160.         Intersection_Found = TRUE;
  161.       }
  162.     }
  163.   }
  164.  
  165.   return(Intersection_Found);
  166. }
  167.  
  168.  
  169.  
  170. /*****************************************************************************
  171. *
  172. * FUNCTION
  173. *
  174. *   Intersect_Quadric
  175. *
  176. * INPUT
  177. *   
  178. * OUTPUT
  179. *   
  180. * RETURNS
  181. *   
  182. * AUTHOR
  183. *
  184. *   POV-Ray Team
  185. *   
  186. * DESCRIPTION
  187. *
  188. *   -
  189. *
  190. * CHANGES
  191. *
  192. *   -
  193. *
  194. ******************************************************************************/
  195.  
  196. static int Intersect_Quadric(Ray, Quadric, Depth1, Depth2)
  197. RAY *Ray;
  198. QUADRIC *Quadric;
  199. DBL *Depth1, *Depth2;
  200. {
  201.   register DBL a, b, c, d;
  202.  
  203.   Increase_Counter(stats[Ray_Quadric_Tests]);
  204.  
  205.   a = Xd * (QA * Xd + QB * Yd + QC * Zd) +
  206.                 Yd * (QE * Yd + QF * Zd) +
  207.                 QH * Zd * Zd;
  208.  
  209.   b = Xd * (QA * Xo + 0.5 * (QB * Yo + QC * Zo + QD)) +
  210.                 Yd * (QE * Yo + 0.5 * (QB * Xo + QF * Zo + QG)) +
  211.                 Zd * (QH * Zo + 0.5 * (QC * Xo + QF * Yo + QI));
  212.  
  213.   c = Xo * (QA * Xo + QB * Yo + QC * Zo + QD) +
  214.                   Yo * (QE * Yo + QF * Zo + QG) +
  215.       Zo * (QH * Zo + QI) +
  216.       QJ;
  217.  
  218.   if (a != 0.0)
  219.   {
  220.     /* The equation is quadratic - find its roots */
  221.  
  222.     d = Sqr(b) - a * c;
  223.  
  224.     if (d <= 0.0)
  225.     {
  226.       return(FALSE);
  227.     }
  228.  
  229.     d = sqrt (d);
  230.  
  231.     *Depth1 = (-b + d) / (a);
  232.     *Depth2 = (-b - d) / (a);
  233.   }
  234.   else
  235.   {
  236.     /* There are no quadratic terms. Solve the linear equation instead. */
  237.  
  238.     if (b == 0.0)
  239.     {
  240.       return(FALSE);
  241.     }
  242.  
  243.     *Depth1 = - 0.5 * c / b;
  244.     *Depth2 = Max_Distance;
  245.   }
  246.  
  247.   Increase_Counter(stats[Ray_Quadric_Tests_Succeeded]);
  248.  
  249.   return(TRUE);
  250. }
  251.  
  252.  
  253.  
  254. /*****************************************************************************
  255. *
  256. * FUNCTION
  257. *
  258. *   Inside_Quadric
  259. *
  260. * INPUT
  261. *   
  262. * OUTPUT
  263. *   
  264. * RETURNS
  265. *   
  266. * AUTHOR
  267. *
  268. *   POV-Ray Team
  269. *   
  270. * DESCRIPTION
  271. *
  272. *   -
  273. *
  274. * CHANGES
  275. *
  276. *   -
  277. *
  278. ******************************************************************************/
  279.  
  280. static int Inside_Quadric(IPoint, Object)
  281. VECTOR IPoint;
  282. OBJECT *Object;
  283. {
  284.   QUADRIC *Quadric = (QUADRIC *)Object;
  285.  
  286.   /* This is faster and shorter. [DB 7/94] */
  287.  
  288.   return((IPoint[X] * (QA * IPoint[X] + QB * IPoint[Y] + QD) +
  289.           IPoint[Y] * (QE * IPoint[Y] + QF * IPoint[Z] + QG) +
  290.           IPoint[Z] * (QH * IPoint[Z] + QC * IPoint[X] + QI) + QJ) <= 0.0);
  291. }
  292.  
  293.  
  294.  
  295. /*****************************************************************************
  296. *
  297. * FUNCTION
  298. *
  299. *   Quadric_Normal
  300. *
  301. * INPUT
  302. *   
  303. * OUTPUT
  304. *   
  305. * RETURNS
  306. *   
  307. * AUTHOR
  308. *
  309. *   POV-Ray Team
  310. *   
  311. * DESCRIPTION
  312. *
  313. *   -
  314. *
  315. * CHANGES
  316. *
  317. *   -
  318. *
  319. ******************************************************************************/
  320.  
  321. static void Quadric_Normal(Result, Object, Inter)
  322. OBJECT *Object;
  323. VECTOR Result;
  324. INTERSECTION *Inter;
  325. {
  326.   QUADRIC *Quadric = (QUADRIC *) Object;
  327.   DBL Len;
  328.  
  329.   /* This is faster and shorter. [DB 7/94] */
  330.  
  331.   Result[X] = 2.0 * QA * Inter->IPoint[X] +
  332.                     QB * Inter->IPoint[Y] +
  333.                     QC * Inter->IPoint[Z] +
  334.                     QD;
  335.  
  336.   Result[Y] =       QB * Inter->IPoint[X] +
  337.               2.0 * QE * Inter->IPoint[Y] +
  338.                     QF * Inter->IPoint[Z] +
  339.                     QG;
  340.  
  341.   Result[Z] =       QC * Inter->IPoint[X] +
  342.                     QF * Inter->IPoint[Y] +
  343.               2.0 * QH * Inter->IPoint[Z] +
  344.                     QI;
  345.  
  346.   VLength(Len, Result);
  347.  
  348.   if (Len == 0.0)
  349.   {
  350.     /* The normal is not defined at this point of the surface. */
  351.     /* Set it to any arbitrary direction. */
  352.  
  353.     Make_Vector(Result, 1.0, 0.0, 0.0);
  354.   }
  355.   else
  356.   {
  357.     VInverseScaleEq(Result, Len);
  358.   }
  359. }
  360.  
  361.  
  362.  
  363. /*****************************************************************************
  364. *
  365. * FUNCTION
  366. *
  367. *   Transform_Quadric
  368. *
  369. * INPUT
  370. *   
  371. * OUTPUT
  372. *   
  373. * RETURNS
  374. *   
  375. * AUTHOR
  376. *
  377. *   POV-Ray Team
  378. *   
  379. * DESCRIPTION
  380. *
  381. *   -
  382. *
  383. * CHANGES
  384. *
  385. *   -
  386. *
  387. ******************************************************************************/
  388.  
  389. static void Transform_Quadric(Object, Trans)
  390. OBJECT *Object;
  391. TRANSFORM *Trans;
  392. {
  393.   QUADRIC *Quadric=(QUADRIC *)Object;
  394.   MATRIX Quadric_Matrix, Transform_Transposed;
  395.  
  396.   Quadric_To_Matrix (Quadric, Quadric_Matrix);
  397.  
  398.   MTimes (Quadric_Matrix, Trans->inverse, Quadric_Matrix);
  399.   MTranspose (Transform_Transposed, Trans->inverse);
  400.   MTimes (Quadric_Matrix, Quadric_Matrix, Transform_Transposed);
  401.  
  402.   Matrix_To_Quadric (Quadric_Matrix, Quadric);
  403.  
  404.   Recompute_BBox(&Object->BBox, Trans);
  405. }
  406.  
  407.  
  408.  
  409. /*****************************************************************************
  410. *
  411. * FUNCTION
  412. *
  413. *   Quadric_To_Matrix
  414. *
  415. * INPUT
  416. *   
  417. * OUTPUT
  418. *   
  419. * RETURNS
  420. *   
  421. * AUTHOR
  422. *
  423. *   POV-Ray Team
  424. *   
  425. * DESCRIPTION
  426. *
  427. *   -
  428. *
  429. * CHANGES
  430. *
  431. *   -
  432. *
  433. ******************************************************************************/
  434.  
  435. static void Quadric_To_Matrix(Quadric, Matrix)
  436. QUADRIC *Quadric;
  437. MATRIX Matrix;
  438. {
  439.   MZero (Matrix);
  440.  
  441.   Matrix[0][0] = Quadric->Square_Terms[X];
  442.   Matrix[1][1] = Quadric->Square_Terms[Y];
  443.   Matrix[2][2] = Quadric->Square_Terms[Z];
  444.   Matrix[0][1] = Quadric->Mixed_Terms[X];
  445.   Matrix[0][2] = Quadric->Mixed_Terms[Y];
  446.   Matrix[0][3] = Quadric->Terms[X];
  447.   Matrix[1][2] = Quadric->Mixed_Terms[Z];
  448.   Matrix[1][3] = Quadric->Terms[Y];
  449.   Matrix[2][3] = Quadric->Terms[Z];
  450.   Matrix[3][3] = Quadric->Constant;
  451. }
  452.  
  453.  
  454.  
  455. /*****************************************************************************
  456. *
  457. * FUNCTION
  458. *
  459. *   Matrix_To_Quadric
  460. *
  461. * INPUT
  462. *   
  463. * OUTPUT
  464. *
  465. * RETURNS
  466. *   
  467. * AUTHOR
  468. *
  469. *   POV-Ray Team
  470. *   
  471. * DESCRIPTION
  472. *
  473. *   -
  474. *
  475. * CHANGES
  476. *
  477. *   -
  478. *
  479. ******************************************************************************/
  480.  
  481. static void Matrix_To_Quadric(Matrix, Quadric)
  482. MATRIX Matrix;
  483. QUADRIC *Quadric;
  484. {
  485.   Quadric->Square_Terms[X] = Matrix[0][0];
  486.   Quadric->Square_Terms[Y] = Matrix[1][1];
  487.   Quadric->Square_Terms[Z] = Matrix[2][2];
  488.  
  489.   Quadric->Mixed_Terms[X] = Matrix[0][1] + Matrix[1][0];
  490.   Quadric->Mixed_Terms[Y] = Matrix[0][2] + Matrix[2][0];
  491.   Quadric->Mixed_Terms[Z] = Matrix[1][2] + Matrix[2][1];
  492.  
  493.   Quadric->Terms[X] = Matrix[0][3] + Matrix[3][0];
  494.   Quadric->Terms[Y] = Matrix[1][3] + Matrix[3][1];
  495.   Quadric->Terms[Z] = Matrix[2][3] + Matrix[3][2];
  496.  
  497.   Quadric->Constant = Matrix[3][3];
  498. }
  499.  
  500.  
  501.  
  502. /*****************************************************************************
  503. *
  504. * FUNCTION
  505. *
  506. *   Translate_Quadric
  507. *
  508. * INPUT
  509. *   
  510. * OUTPUT
  511. *   
  512. * RETURNS
  513. *   
  514. * AUTHOR
  515. *
  516. *   POV-Ray Team
  517. *   
  518. * DESCRIPTION
  519. *
  520. *   -
  521. *
  522. * CHANGES
  523. *
  524. *   -
  525. *
  526. ******************************************************************************/
  527.  
  528. static void Translate_Quadric(Object, Vector, Trans)
  529. OBJECT *Object;
  530. VECTOR Vector;
  531. TRANSFORM *Trans;
  532. {
  533.   Transform_Quadric (Object, Trans);
  534. }
  535.  
  536.  
  537.  
  538. /*****************************************************************************
  539. *
  540. * FUNCTION
  541. *
  542. *   Rotate_Quadric
  543. *
  544. * INPUT
  545. *   
  546. * OUTPUT
  547. *   
  548. * RETURNS
  549. *   
  550. * AUTHOR
  551. *
  552. *   POV-Ray Team
  553. *   
  554. * DESCRIPTION
  555. *
  556. *   -
  557. *
  558. * CHANGES
  559. *
  560. *   -
  561. *
  562. ******************************************************************************/
  563.  
  564. static void Rotate_Quadric(Object, Vector, Trans)
  565. OBJECT *Object;
  566. VECTOR Vector;
  567. TRANSFORM *Trans;
  568. {
  569.   Transform_Quadric (Object, Trans);
  570. }
  571.  
  572.  
  573.  
  574. /*****************************************************************************
  575. *
  576. * FUNCTION
  577. *
  578. *   Scale_Quadric
  579. *
  580. * INPUT
  581. *   
  582. * OUTPUT
  583. *   
  584. * RETURNS
  585. *   
  586. * AUTHOR
  587. *
  588. *   POV-Ray Team
  589. *   
  590. * DESCRIPTION
  591. *
  592. *   -
  593. *
  594. * CHANGES
  595. *
  596. *   -
  597. *
  598. ******************************************************************************/
  599.  
  600. static void Scale_Quadric(Object, Vector, Trans)
  601. OBJECT *Object;
  602. VECTOR Vector;
  603. TRANSFORM *Trans;
  604. {
  605.   Transform_Quadric (Object, Trans);
  606. }
  607.  
  608.  
  609.  
  610. /*****************************************************************************
  611. *
  612. * FUNCTION
  613. *
  614. *   Invert_Quadric
  615. *
  616. * INPUT
  617. *   
  618. * OUTPUT
  619. *   
  620. * RETURNS
  621. *   
  622. * AUTHOR
  623. *
  624. *   POV-Ray Team
  625. *   
  626. * DESCRIPTION
  627. *
  628. *   -
  629. *
  630. * CHANGES
  631. *
  632. *   -
  633. *
  634. ******************************************************************************/
  635.  
  636. static void Invert_Quadric(Object)
  637. OBJECT *Object;
  638. {
  639.   QUADRIC *Quadric = (QUADRIC *) Object;
  640.  
  641.   VScaleEq(Quadric->Square_Terms, -1.0);
  642.   VScaleEq(Quadric->Mixed_Terms, -1.0);
  643.   VScaleEq(Quadric->Terms, -1.0);
  644.  
  645.   Quadric->Constant *= -1.0;
  646.  
  647.   Invert_Flag(Object, INVERTED_FLAG);
  648. }
  649.  
  650.  
  651.  
  652. /*****************************************************************************
  653. *
  654. * FUNCTION
  655. *
  656. *   Create_Quadric
  657. *
  658. * INPUT
  659. *   
  660. * OUTPUT
  661. *   
  662. * RETURNS
  663. *   
  664. * AUTHOR
  665. *
  666. *   POV-Ray Team
  667. *   
  668. * DESCRIPTION
  669. *
  670. *   -
  671. *
  672. * CHANGES
  673. *
  674. *   -
  675. *
  676. ******************************************************************************/
  677.  
  678. QUADRIC *Create_Quadric()
  679. {
  680.   QUADRIC *New;
  681.  
  682.   New = (QUADRIC *)POV_MALLOC(sizeof (QUADRIC), "quadric");
  683.  
  684.   INIT_OBJECT_FIELDS(New, QUADRIC_OBJECT, &Quadric_Methods)
  685.  
  686.   Make_Vector (New->Square_Terms, 1.0, 1.0, 1.0);
  687.   Make_Vector (New->Mixed_Terms, 0.0, 0.0, 0.0);
  688.   Make_Vector (New->Terms, 0.0, 0.0, 0.0);
  689.   New->Constant = 1.0;
  690.  
  691.   return(New);
  692. }
  693.  
  694.  
  695.  
  696. /*****************************************************************************
  697. *
  698. * FUNCTION
  699. *
  700. *   Copy_Quadric
  701. *
  702. * INPUT
  703. *   
  704. * OUTPUT
  705. *   
  706. * RETURNS
  707. *   
  708. * AUTHOR
  709. *
  710. *   POV-Ray Team
  711. *   
  712. * DESCRIPTION
  713. *
  714. *   -
  715. *
  716. * CHANGES
  717. *
  718. *   -
  719. *
  720. ******************************************************************************/
  721.  
  722. static void *Copy_Quadric(Object)
  723. OBJECT *Object;
  724. {
  725.   QUADRIC *New;
  726.  
  727.   New = Create_Quadric();
  728.  
  729.   *New = *((QUADRIC *)Object);
  730.  
  731.   return(New);
  732. }
  733.  
  734.  
  735.  
  736. /*****************************************************************************
  737. *
  738. * FUNCTION
  739. *
  740. *   Destroy_Quadric
  741. *
  742. * INPUT
  743. *   
  744. * OUTPUT
  745. *   
  746. * RETURNS
  747. *   
  748. * AUTHOR
  749. *
  750. *   POV-Ray Team
  751. *   
  752. * DESCRIPTION
  753. *
  754. *   -
  755. *
  756. * CHANGES
  757. *
  758. *   -
  759. *
  760. ******************************************************************************/
  761.  
  762. static void Destroy_Quadric(Object)
  763. OBJECT *Object;
  764. {
  765.   POV_FREE (Object);
  766. }
  767.  
  768.  
  769.  
  770. /*****************************************************************************
  771. *
  772. * FUNCTION
  773. *
  774. *   Compute_Quadric_BBox
  775. *
  776. * INPUT
  777. *
  778. *   Quadric - Qaudric object
  779. *   
  780. * OUTPUT
  781. *
  782. *   Quadric
  783. *   
  784. * RETURNS
  785. *   
  786. * AUTHOR
  787. *
  788. *   Dieter Bayer
  789. *   
  790. * DESCRIPTION
  791. *
  792. *   Compute a bounding box around a quadric.
  793. *
  794. *   This function calculates the bounding box of a quadric given in
  795. *   its normal form, i.e. f(x,y,z) = A*x^2 + E*y^2 + H*z^2 + J = 0.
  796. *
  797. *   NOTE: Translated quadrics can also be bounded by this function.
  798. *
  799. *         Supported: cones, cylinders, ellipsoids, hyperboloids.
  800. *
  801. * CHANGES
  802. *
  803. *   May 1994 : Creation.
  804. *
  805. *   Sep 1994 : Added support of hyperpoloids. Improved bounding of
  806. *              quadrics used in CSG intersections. [DB]
  807. *
  808. ******************************************************************************/
  809.  
  810. void Compute_Quadric_BBox(Quadric, ClipMin, ClipMax)
  811. QUADRIC *Quadric;
  812. VECTOR ClipMin, ClipMax;
  813. {
  814.   DBL A, B, C, D, E, F, G, H, I, J;
  815.   DBL a, b, c, d;
  816.   DBL rx, ry, rz, rx1, rx2, ry1, ry2, rz1, rz2, x, y, z;
  817.   DBL New_Volume, Old_Volume;
  818.   VECTOR Min, Max, TmpMin, TmpMax, NewMin, NewMax, T1;
  819.   BBOX Old;
  820.   OBJECT *Sib;
  821.  
  822.   /*
  823.    * Check for 'normal' form. If the quadric isn't in it's normal form
  824.    * we can't do anything (we could, but that would be to tedious!
  825.    * Diagonalising the quadric's 4x4 matrix, i.e. finding its eigenvalues
  826.    * and eigenvectors -> solving a 4th order polynom).
  827.    */
  828.  
  829.   /* Get quadrics coefficients. */
  830.  
  831.   A = Quadric->Square_Terms[X];
  832.   E = Quadric->Square_Terms[Y];
  833.   H = Quadric->Square_Terms[Z];
  834.   B = Quadric->Mixed_Terms[X] / 2.0;
  835.   C = Quadric->Mixed_Terms[Y] / 2.0;
  836.   F = Quadric->Mixed_Terms[Z] / 2.0;
  837.   D = Quadric->Terms[X] / 2.0;
  838.   G = Quadric->Terms[Y] / 2.0;
  839.   I = Quadric->Terms[Z] / 2.0;
  840.   J = Quadric->Constant;
  841.  
  842.   /* Set small values to 0. */
  843.  
  844.   if (fabs(A) < EPSILON) A = 0.0;
  845.   if (fabs(B) < EPSILON) B = 0.0;
  846.   if (fabs(C) < EPSILON) C = 0.0;
  847.   if (fabs(D) < EPSILON) D = 0.0;
  848.   if (fabs(E) < EPSILON) E = 0.0;
  849.   if (fabs(F) < EPSILON) F = 0.0;
  850.   if (fabs(G) < EPSILON) G = 0.0;
  851.   if (fabs(H) < EPSILON) H = 0.0;
  852.   if (fabs(I) < EPSILON) I = 0.0;
  853.   if (fabs(J) < EPSILON) J = 0.0;
  854.  
  855.   /* Non-zero mixed terms --> return */
  856.  
  857.   if ((B != 0.0) || (C != 0.0) || (F != 0.0))
  858.   {
  859.     return;
  860.   }
  861.  
  862.   /* Non-zero linear terms --> get translation vector */
  863.  
  864.   if ((D != 0.0) || (G != 0.0) || (I != 0.0))
  865.   {
  866.     if (A != 0.0)
  867.     {
  868.       T1[X] = -D / A;
  869.     }
  870.     else
  871.     {
  872.       if (D != 0.0)
  873.       {
  874.        T1[X] = J / (2.0 * D);
  875.       }
  876.       else
  877.       {
  878.         T1[X] = 0.0;
  879.       }
  880.     }
  881.  
  882.     if (E != 0.0)
  883.     {
  884.       T1[Y] = -G / E;
  885.     }
  886.     else
  887.     {
  888.       if (G != 0.0)
  889.       {
  890.         T1[Y] = J / (2.0 * G);
  891.       }
  892.       else
  893.       {
  894.         T1[Y] = 0.0;
  895.       }
  896.     }
  897.  
  898.     if (H != 0.0)
  899.     {
  900.       T1[Z] = -I / H;
  901.     }
  902.     else
  903.     {
  904.       if (I != 0.0)
  905.       {
  906.         T1[Z] = J / (2.0 * I);
  907.       }
  908.       else
  909.       {
  910.         T1[Z] = 0.0;
  911.       }
  912.     }
  913.  
  914.     /* Recalculate coefficients. */
  915.  
  916.     D += A * T1[X];
  917.     G += E * T1[Y];
  918.     I += H * T1[Z];
  919.     J -= T1[X]*(A*T1[X] + 2.0*D) + T1[Y]*(E*T1[Y] + 2.0*G) + T1[Z]*(H*T1[Z] + 2.0*I);
  920.   }
  921.   else
  922.   {
  923.     Make_Vector(T1, 0.0, 0.0, 0.0);
  924.   }
  925.  
  926.   /* Get old bounding box. */
  927.  
  928.   Old = Quadric->BBox;
  929.  
  930.   /* Init new bounding box. */
  931.  
  932.   NewMin[X] = NewMin[Y] = NewMin[Z] = -BOUND_HUGE/2;
  933.   NewMax[X] = NewMax[Y] = NewMax[Z] =  BOUND_HUGE/2;
  934.  
  935.   /* Get the bounding box of the clipping object. */
  936.  
  937.   if (Quadric->Clip != NULL)
  938.   {
  939.     Min[X] = Min[Y] = Min[Z] = -BOUND_HUGE;
  940.     Max[X] = Max[Y] = Max[Z] =  BOUND_HUGE;
  941.  
  942.     /* Intersect the members bounding boxes. */
  943.  
  944.     for (Sib = Quadric->Clip; Sib != NULL; Sib = Sib->Sibling)
  945.     {
  946.       if (!Test_Flag(Sib, INVERTED_FLAG))
  947.       {
  948.         if (Sib->Methods == &Plane_Methods)
  949.         {
  950.           Compute_Plane_Min_Max((PLANE *)Sib, TmpMin, TmpMax);
  951.         }
  952.         else
  953.         {
  954.           Make_min_max_from_BBox(TmpMin, TmpMax, Sib->BBox);
  955.         }
  956.  
  957.         Min[X] = max(Min[X], TmpMin[X]);
  958.         Min[Y] = max(Min[Y], TmpMin[Y]);
  959.         Min[Z] = max(Min[Z], TmpMin[Z]);
  960.         Max[X] = min(Max[X], TmpMax[X]);
  961.         Max[Y] = min(Max[Y], TmpMax[Y]);
  962.         Max[Z] = min(Max[Z], TmpMax[Z]);
  963.       }
  964.     }
  965.  
  966.     Assign_Vector(ClipMin, Min);
  967.     Assign_Vector(ClipMax, Max);
  968.   }
  969.  
  970.   /* Translate clipping box. */
  971.  
  972.   VSubEq(ClipMin, T1);
  973.   VSubEq(ClipMax, T1);
  974.  
  975.   /* We want A to be non-negative. */
  976.  
  977.   if (A < 0.0)
  978.   {
  979.     A = -A;
  980.     D = -D;
  981.     E = -E;
  982.     G = -G;
  983.     H = -H;
  984.     I = -I;
  985.     J = -J;
  986.   }
  987.  
  988.   /*
  989.    *
  990.    * Check for ellipsoid.
  991.    *
  992.    *    x*x     y*y     z*z
  993.    *   ----- + ----- + ----- - 1 = 0
  994.    *    a*a     b*b     c*c
  995.    *
  996.    */
  997.  
  998.   if ((A > 0.0) && (E > 0.0) && (H > 0.0) && (J < 0.0))
  999.   {
  1000.     a = sqrt(-J/A);
  1001.     b = sqrt(-J/E);
  1002.     c = sqrt(-J/H);
  1003.  
  1004.     NewMin[X] = -a;
  1005.     NewMin[Y] = -b;
  1006.     NewMin[Z] = -c;
  1007.     NewMax[X] = a;
  1008.     NewMax[Y] = b;
  1009.     NewMax[Z] = c;
  1010.   }
  1011.  
  1012.   /*
  1013.    *
  1014.    * Check for cylinder (x-axis).
  1015.    *
  1016.    *    y*y     z*z
  1017.    *   ----- + ----- - 1 = 0
  1018.    *    b*b     c*c
  1019.    *
  1020.    */
  1021.  
  1022.   if ((A == 0.0) && (E > 0.0) && (H > 0.0) && (J < 0.0))
  1023.   {
  1024.     b = sqrt(-J/E);
  1025.     c = sqrt(-J/H);
  1026.  
  1027.     NewMin[Y] = -b;
  1028.     NewMin[Z] = -c;
  1029.     NewMax[Y] = b;
  1030.     NewMax[Z] = c;
  1031.   }
  1032.  
  1033.   /*
  1034.    *
  1035.    * Check for cylinder (y-axis).
  1036.    *
  1037.    *    x*x     z*z
  1038.    *   ----- + ----- - 1 = 0
  1039.    *    a*a     c*c
  1040.    *
  1041.    */
  1042.  
  1043.   if ((A > 0.0) && (E == 0.0) && (H > 0.0) && (J < 0.0))
  1044.   {
  1045.     a = sqrt(-J/A);
  1046.     c = sqrt(-J/H);
  1047.  
  1048.     NewMin[X] = -a;
  1049.     NewMin[Z] = -c;
  1050.     NewMax[X] = a;
  1051.     NewMax[Z] = c;
  1052.   }
  1053.  
  1054.   /*
  1055.    *
  1056.    * Check for cylinder (z-axis).
  1057.    * 
  1058.    *    x*x     y*y
  1059.    *   ----- + ----- - 1 = 0
  1060.    *    a*a     b*b
  1061.    *
  1062.    */
  1063.  
  1064.   if ((A > 0.0) && (E > 0.0) && (H == 0.0) && (J < 0.0))
  1065.   {
  1066.     a = sqrt(-J/A);
  1067.     b = sqrt(-J/E);
  1068.  
  1069.     NewMin[X] = -a;
  1070.     NewMin[Y] = -b;
  1071.     NewMax[X] = a;
  1072.     NewMax[Y] = b;
  1073.   }
  1074.  
  1075.   /*
  1076.    *
  1077.    * Check for cone (x-axis).
  1078.    *
  1079.    *    x*x     y*y     z*z
  1080.    *   ----- - ----- - ----- = 0
  1081.    *    a*a     b*b     c*c
  1082.    *
  1083.    */
  1084.  
  1085.   if ((A > 0.0) && (E < 0.0) && (H < 0.0) && (J == 0.0))
  1086.   {
  1087.     a = sqrt(1.0/A);
  1088.     b = sqrt(-1.0/E);
  1089.     c = sqrt(-1.0/H);
  1090.  
  1091.     /* Get radii for lower x value. */
  1092.  
  1093.     x = ClipMin[X];
  1094.  
  1095.     ry1 = fabs(x * b / a);
  1096.     rz1 = fabs(x * c / a);
  1097.  
  1098.     /* Get radii for upper x value. */
  1099.  
  1100.     x = ClipMax[X];
  1101.  
  1102.     ry2 = fabs(x * b / a);
  1103.     rz2 = fabs(x * c / a);
  1104.  
  1105.     ry = max(ry1, ry2);
  1106.     rz = max(rz1, rz2);
  1107.  
  1108.     NewMin[Y] = -ry;
  1109.     NewMin[Z] = -rz;
  1110.     NewMax[Y] = ry;
  1111.     NewMax[Z] = rz;
  1112.   }
  1113.  
  1114.   /*
  1115.    *
  1116.    *  Check for cone (y-axis).
  1117.    *
  1118.    *    x*x     y*y     z*z
  1119.    *   ----- - ----- + ----- = 0
  1120.    *    a*a     b*b     c*c
  1121.    *
  1122.    */
  1123.  
  1124.   if ((A > 0.0) && (E < 0.0) && (H > 0.0) && (J == 0.0))
  1125.   {
  1126.     a = sqrt(1.0/A);
  1127.     b = sqrt(-1.0/E);
  1128.     c = sqrt(1.0/H);
  1129.  
  1130.     /* Get radii for lower y value. */
  1131.  
  1132.     y = ClipMin[Y];
  1133.  
  1134.     rx1 = fabs(y * a / b);
  1135.     rz1 = fabs(y * c / b);
  1136.  
  1137.     /* Get radii for upper y value. */
  1138.  
  1139.     y = ClipMax[Y];
  1140.  
  1141.     rx2 = fabs(y * a / b);
  1142.     rz2 = fabs(y * c / b);
  1143.  
  1144.     rx = max(rx1, rx2);
  1145.     rz = max(rz1, rz2);
  1146.  
  1147.     NewMin[X] = -rx;
  1148.     NewMin[Z] = -rz;
  1149.     NewMax[X] = rx;
  1150.     NewMax[Z] = rz;
  1151.   }
  1152.  
  1153.   /*
  1154.    *
  1155.    * Check for cone (z-axis).
  1156.    * 
  1157.    *    x*x     y*y     z*z
  1158.    *   ----- + ----- - ----- = 0
  1159.    *    a*a     b*b     c*c
  1160.    *
  1161.    */
  1162.  
  1163.   if ((A > 0.0) && (E > 0.0) && (H < 0.0) && (J == 0.0))
  1164.   {
  1165.     a = sqrt(1.0/A);
  1166.     b = sqrt(1.0/E);
  1167.     c = sqrt(-1.0/H);
  1168.  
  1169.     /* Get radii for lower z value. */
  1170.  
  1171.     z = ClipMin[Z];
  1172.  
  1173.     rx1 = fabs(z * a / c);
  1174.     ry1 = fabs(z * b / c);
  1175.  
  1176.     /* Get radii for upper z value. */
  1177.  
  1178.     z = ClipMax[Z];
  1179.  
  1180.     rx2 = fabs(z * a / c);
  1181.     ry2 = fabs(z * b / c);
  1182.  
  1183.     rx = max(rx1, rx2);
  1184.     ry = max(ry1, ry2);
  1185.  
  1186.     NewMin[X] = -rx;
  1187.     NewMin[Y] = -ry;
  1188.     NewMax[X] = rx;
  1189.     NewMax[Y] = ry;
  1190.   }
  1191.  
  1192.   /*
  1193.    *
  1194.    * Check for hyperboloid (x-axis).
  1195.    *
  1196.    *    x*x     y*y     z*z
  1197.    *   ----- - ----- - ----- + 1 = 0
  1198.    *    a*a     b*b     c*c
  1199.    *
  1200.    */
  1201.  
  1202.   if ((A > 0.0) && (E < 0.0) && (H < 0.0) && (J > 0.0))
  1203.   {
  1204.     /* Get radii for lower x value. */
  1205.  
  1206.     x = ClipMin[X];
  1207.  
  1208.     d = 1.0 + A * Sqr(x);
  1209.  
  1210.     ry1 = sqrt(-d / E);
  1211.     rz1 = sqrt(-d / H);
  1212.  
  1213.     /* Get radii for upper x value. */
  1214.  
  1215.     x = ClipMax[X];
  1216.  
  1217.     d = 1.0 + A * Sqr(x);
  1218.  
  1219.     ry2 = sqrt(-d / E);
  1220.     rz2 = sqrt(-d / H);
  1221.  
  1222.     ry = max(ry1, ry2);
  1223.     rz = max(rz1, rz2);
  1224.  
  1225.     NewMin[Y] = -ry;
  1226.     NewMin[Z] = -rz;
  1227.     NewMax[Y] = ry;
  1228.     NewMax[Z] = rz;
  1229.   }
  1230.  
  1231.   /*
  1232.    *
  1233.    * Check for hyperboloid (y-axis).
  1234.    * 
  1235.    *    x*x     y*y     z*z
  1236.    *   ----- - ----- + ----- - 1 = 0
  1237.    *    a*a     b*b     c*c
  1238.    *
  1239.    */
  1240.  
  1241.   if ((A > 0.0) && (E < 0.0) && (H > 0.0) && (J < 0.0))
  1242.   {
  1243.     /* Get radii for lower y value. */
  1244.  
  1245.     y = ClipMin[Y];
  1246.  
  1247.     d = 1.0 - E * Sqr(y);
  1248.  
  1249.     rx1 = sqrt(d / A);
  1250.     rz1 = sqrt(d / H);
  1251.  
  1252.     /* Get radii for upper y value. */
  1253.  
  1254.     y = ClipMax[Y];
  1255.  
  1256.     d = 1.0 - E * Sqr(y);
  1257.  
  1258.     rx2 = sqrt(d / A);
  1259.     rz2 = sqrt(d / H);
  1260.  
  1261.     rx = max(rx1, rx2);
  1262.     rz = max(rz1, rz2);
  1263.  
  1264.     NewMin[X] = -rx;
  1265.     NewMin[Z] = -rz;
  1266.     NewMax[X] = rx;
  1267.     NewMax[Z] = rz;
  1268.   }
  1269.  
  1270.   /*
  1271.    *
  1272.    * Check for hyperboloid (z-axis).
  1273.    *
  1274.    *    x*x     y*y     z*z
  1275.    *   ----- + ----- - ----- - 1 = 0
  1276.    *    a*a     b*b     c*c
  1277.    *
  1278.    */
  1279.  
  1280.   if ((A > 0.0) && (E > 0.0) && (H < 0.0) && (J < 0.0))
  1281.   {
  1282.     /* Get radii for lower z value. */
  1283.  
  1284.     z = ClipMin[Z];
  1285.  
  1286.     d = 1.0 - H * Sqr(z);
  1287.  
  1288.     rx1 = sqrt(d / A);
  1289.     ry1 = sqrt(d / E);
  1290.  
  1291.     /* Get radii for upper z value. */
  1292.  
  1293.     z = ClipMax[Z];
  1294.  
  1295.     d = 1.0 - H * Sqr(z);
  1296.  
  1297.     rx2 = sqrt(d / A);
  1298.     ry2 = sqrt(d / E);
  1299.  
  1300.     rx = max(rx1, rx2);
  1301.     ry = max(ry1, ry2);
  1302.  
  1303.     NewMin[X] = -rx;
  1304.     NewMin[Y] = -ry;
  1305.     NewMax[X] = rx;
  1306.     NewMax[Y] = ry;
  1307.   }
  1308.  
  1309.   /*
  1310.    *
  1311.    * Check for paraboloid (x-axis).
  1312.    *
  1313.    *        y*y     z*z
  1314.    *   x - ----- - ----- = 0
  1315.    *        b*b     c*c
  1316.    *
  1317.    */
  1318.  
  1319.   if ((A == 0.0) && (D != 0.0) && (E != 0.0) && (H != 0.0) && (J == 0.0))
  1320.   {
  1321.     /* Get radii for lower x value. */
  1322.  
  1323.     x = ClipMin[X];
  1324.  
  1325.     ry1 = sqrt(fabs(2.0 * D * x / E));
  1326.     rz1 = sqrt(fabs(2.0 * D * x / H));
  1327.  
  1328.     /* Get radii for upper x value. */
  1329.  
  1330.     x = ClipMax[X];
  1331.  
  1332.     ry2 = sqrt(fabs(2.0 * D * x / E));
  1333.     rz2 = sqrt(fabs(2.0 * D * x / H));
  1334.  
  1335.     ry = max(ry1, ry2);
  1336.     rz = max(rz1, rz2);
  1337.  
  1338.     NewMin[Y] = -ry;
  1339.     NewMin[Z] = -rz;
  1340.     NewMax[Y] = ry;
  1341.     NewMax[Z] = rz;
  1342.   }
  1343.  
  1344.   /*
  1345.    *
  1346.    * Check for paraboloid (y-axis).
  1347.    *
  1348.    *        x*x     z*z
  1349.    *   y - ----- - ----- = 0
  1350.    *        a*a     c*c
  1351.    *
  1352.    */
  1353.  
  1354.   if ((E == 0.0) && (G != 0.0) && (A != 0.0) && (H != 0.0) && (J == 0.0))
  1355.   {
  1356.     /* Get radii for lower y-value. */
  1357.  
  1358.     y = ClipMin[Y];
  1359.  
  1360.     rx1 = sqrt(fabs(2.0 * G * y / A));
  1361.     rz1 = sqrt(fabs(2.0 * G * y / H));
  1362.  
  1363.     /* Get radii for upper y value. */
  1364.  
  1365.     y = ClipMax[Y];
  1366.  
  1367.     rx2 = sqrt(fabs(2.0 * G * y / A));
  1368.     rz2 = sqrt(fabs(2.0 * G * y / H));
  1369.  
  1370.     rx = max(rx1, rx2);
  1371.     rz = max(rz1, rz2);
  1372.  
  1373.     NewMin[X] = -rx;
  1374.     NewMin[Z] = -rz;
  1375.     NewMax[X] = rx;
  1376.     NewMax[Z] = rz;
  1377.   }
  1378.  
  1379.   /*
  1380.    *
  1381.    * Check for paraboloid (z-axis).
  1382.    *
  1383.    *        x*x     y*y
  1384.    *   z - ----- - ----- = 0
  1385.    *        a*a     b*b
  1386.    *
  1387.    */
  1388.  
  1389.   if ((H == 0.0) && (I != 0.0) && (A != 0.0) && (E != 0.0) && (J == 0.0))
  1390.   {
  1391.     /* Get radii for lower z-value. */
  1392.  
  1393.     z = ClipMin[Z];
  1394.  
  1395.     rx1 = sqrt(fabs(2.0 * I * z / A));
  1396.     ry1 = sqrt(fabs(2.0 * I * z / E));
  1397.  
  1398.     /* Get radii for upper z value. */
  1399.  
  1400.     z = ClipMax[Z];
  1401.  
  1402.     rx2 = sqrt(fabs(2.0 * I * z / A));
  1403.     ry2 = sqrt(fabs(2.0 * I * z / E));
  1404.  
  1405.     rx = max(rx1, rx2);
  1406.     ry = max(ry1, ry2);
  1407.  
  1408.     NewMin[X] = -rx;
  1409.     NewMin[Y] = -ry;
  1410.     NewMax[X] = rx;
  1411.     NewMax[Y] = ry;
  1412.   }
  1413.  
  1414.   /* Intersect clipping object's and quadric's bounding boxes */
  1415.  
  1416.   NewMin[X] = max(NewMin[X], ClipMin[X]);
  1417.   NewMin[Y] = max(NewMin[Y], ClipMin[Y]);
  1418.   NewMin[Z] = max(NewMin[Z], ClipMin[Z]);
  1419.  
  1420.   NewMax[X] = min(NewMax[X], ClipMax[X]);
  1421.   NewMax[Y] = min(NewMax[Y], ClipMax[Y]);
  1422.   NewMax[Z] = min(NewMax[Z], ClipMax[Z]);
  1423.  
  1424.   /* Use old or new bounding box? */
  1425.  
  1426.   New_Volume = (NewMax[X] - NewMin[X]) * (NewMax[Y] - NewMin[Y]) * (NewMax[Z] - NewMin[Z]);
  1427.  
  1428.   BOUNDS_VOLUME(Old_Volume, Old);
  1429.  
  1430.   if (New_Volume < Old_Volume)
  1431.   {
  1432.     /* Add translation. */
  1433.  
  1434.     VAddEq(NewMin, T1);
  1435.     VAddEq(NewMax, T1);
  1436.  
  1437.     Make_BBox_from_min_max(Quadric->BBox, NewMin, NewMax);
  1438.  
  1439.     /* Beware of bounding boxes to large. */
  1440.  
  1441.     if ((Quadric->BBox.Lengths[X] > CRITICAL_LENGTH) ||
  1442.         (Quadric->BBox.Lengths[Y] > CRITICAL_LENGTH) ||
  1443.         (Quadric->BBox.Lengths[Z] > CRITICAL_LENGTH))
  1444.     {
  1445.       Make_BBox(Quadric->BBox, -BOUND_HUGE/2, -BOUND_HUGE/2, -BOUND_HUGE/2,
  1446.         BOUND_HUGE, BOUND_HUGE, BOUND_HUGE);
  1447.     }
  1448.   }
  1449. }
  1450.  
  1451.  
  1452.  
  1453. /*****************************************************************************
  1454. *
  1455. * FUNCTION
  1456. *
  1457. *   Compute_Plane_Min_Max
  1458. *
  1459. * INPUT
  1460. *
  1461. *   Plane    - Plane
  1462. *   Min, Max - Vectors containing plane's dimensions
  1463. *   
  1464. * OUTPUT
  1465. *
  1466. *   Min, Max
  1467. *   
  1468. * RETURNS
  1469. *   
  1470. * AUTHOR
  1471. *
  1472. *   Dieter Bayer
  1473. *   
  1474. * DESCRIPTION
  1475. *
  1476. *   Compute min/max vectors for planes perpendicular to an axis.
  1477. *
  1478. * CHANGES
  1479. *
  1480. *   May 1994 : Creation.
  1481. *
  1482. ******************************************************************************/
  1483.  
  1484. void Compute_Plane_Min_Max(Plane, Min, Max)
  1485. PLANE *Plane;
  1486. VECTOR Min, Max;
  1487. {
  1488.   DBL d;
  1489.   VECTOR N;
  1490.  
  1491.   Assign_Vector(N, Plane->Normal_Vector);
  1492.  
  1493.   d = -Plane->Distance;
  1494.  
  1495.   Min[X] = Min[Y] = Min[Z] = -BOUND_HUGE/2;
  1496.   Max[X] = Max[Y] = Max[Z] =  BOUND_HUGE/2;
  1497.  
  1498.   /* y-z-plane */
  1499.  
  1500.   if (fabs(1.0 - fabs(N[X])) < EPSILON)
  1501.   {
  1502.     if (N[X] > 0.0)
  1503.     {
  1504.       Max[X] = d;
  1505.     }
  1506.     else
  1507.     {
  1508.       Min[X] = -d;
  1509.     }
  1510.   }
  1511.  
  1512.   /* x-z-plane */
  1513.  
  1514.   if (fabs(1.0 - fabs(N[Y])) < EPSILON)
  1515.   {
  1516.     if (N[Y] > 0.0)
  1517.     {
  1518.       Max[Y] = d;
  1519.     }
  1520.     else
  1521.     {
  1522.       Min[Y] = -d;
  1523.     }
  1524.   }
  1525.  
  1526.   /* x-y-plane */
  1527.  
  1528.   if (fabs(1.0 - fabs(N[Z])) < EPSILON)
  1529.   {
  1530.     if (N[Z] > 0.0)
  1531.     {
  1532.       Max[Z] = d;
  1533.     }
  1534.     else
  1535.     {
  1536.       Min[Z] = -d;
  1537.     }
  1538.   }
  1539. }
  1540.  
  1541.  
  1542.  
  1543.